###maps_create.r
####create all maps, both static and otherwise
#####RdE 11/2012

rm(list = ls())


##load necessary packages
library(maptools)
library(raster)
library(colorspace)
library(rgdal)

cat('loading data...\n\n')

##us data, from web source
us1 <- getData('GADM', country="USA", level=1)
us2 <- getData('GADM', country="USA", level=2)

##dma data, proprietary from ArcGIS
dma.state <-readShapePoly("R_map_files/state_DMA_intersections.shp")
dma.use <-readShapePoly("R_map_files/DMA_with_variation.shp")

##create new identifier to keep track of state/dma combinations
dma.state$dma.state = paste(dma.state$DMA_2012_I,dma.state$US_state_2, sep = '')


##independent variables
###limit romney effort states to states that have DMA's in sample
##romney state data
romney.state = read.csv('R_map_files/romney_effort.csv')
romney.state = romney.state[romney.state$romney_effort>0,]
##order for use below
romney.state = romney.state[order(romney.state$romney_effort),]

##battleground states
bgs = c('Colorado', 'Florida', 'Iowa', 'Nevada', 'New Hampshire', 'North Carolina', 'Ohio', 'Virginia',
        'Wisconsin') 
dma.state$bg = ifelse(dma.state$US_state_3 %in% bgs,1,0)

bg.dmas = c('790','751','561','686','530','682','652','631','717','611','624','725','811','770',
            '506','523','500','524','517','575','570','518','567','544','560','564','515','542',
            '509','597','547','554','536','559','518','569','544','560','531','676','658','702',
            '553','613')

##tiers
tier1 = c('Colorado', 'Iowa', 'Nevada', 'New Hampshire', 'Ohio', 'Virginia')
tier2 = c('Pennsylvania', 'Florida', 'North Carolina')
tier3 = c('Wisconsin', 'Minnesota',  'Oregon', 'Michigan', 'New Mexico')

tier.dmas = c('790','751','561','686','530','682','652','631','717','611','624','725','811','770',
              '506','523','500','524','517','575','570','518','567','544','560','564','515','542',
              '509','597','547','554','536','559','518','569','544','560','531','676','658','702',
              '553','613','676','724','702','613','611','725','813','676','658','553','588','547',
              '690','634','765','514','565','501','504','508','511','536')


##modify spatial data
##add projection to DMA's 
##for dmas
proj4string(dma.state) <- CRS("+proj=merc +a=6378137 +b=6356752.314245179 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +ellps=WGS84 +datum=WGS84 +units=m +nadgrids=@null +wktext +no_defs")
proj4string(dma.use) <- CRS("+proj=merc +a=6378137 +b=6356752.314245179 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0.0 +k=1.0 +ellps=WGS84 +datum=WGS84 +units=m +nadgrids=@null +wktext +no_defs")

##extract Michigan from counties to use below, remove bodies of water
mi2 = us2[us2$NAME_1 == 'Michigan' & us2$TYPE_2!='Water body',]
##remove michigan because of odd shape, use trick below to recombine.  
us1 = us1[us1$NAME_1 != 'Michigan',]
##merge the Michigan counties into a single new state, which makes it easier to work with below
mi2 = unionSpatialPolygons(mi2,ID = mi2$ID_1)

## Specify a geographic extent for the map
## by defining the top-left and bottom-right geographic coordinates
left = -118; top = 50; right = -83; bottom = 22
mapExtent <- rbind(c(left, top), c(right,bottom))

## Specify the required projection using a proj4 string
## Use http://www.spatialreference.org/ to find the required string
## Polyconic for North America
newProj <- CRS("+proj=poly +lat_0=0 +lon_0=-100 +x_0=0 
  		+y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")


## Project the map extent (first need to specify that it is longlat) 
nonproj.coords = SpatialPoints(mapExtent, proj4string=CRS("+proj=longlat"))
mapExtentPr <- spTransform(SpatialPoints(mapExtent, 
				  proj4string=CRS("+proj=longlat")),
				  newProj)

##get these points out of the coordinate system for use in drawing rectangle below
top.proj = bbox(mapExtentPr)['coords.x2','max'] 
left.proj = bbox(mapExtentPr)['coords.x1','min']
bottom.proj = bbox(mapExtentPr)['coords.x2','min'] 
right.proj = bbox(mapExtentPr)['coords.x1','max'] 		  

us1Pr <- spTransform(us1, newProj) 
mi2Pr <- spTransform(mi2, newProj) 
dma.state <- spTransform(dma.state, newProj) 
dma.use <- spTransform(dma.use, newProj) 


##colors for maps
###set color elements for use below
##use this for mapping places with zero romney effort
zero.color = 'gray'  ##color to represent zero below
campaign.chroma = c(32, 50)
campaign.lum = c(80,16)

##values of Romney effort, will be matched with colors
romney.effort = unique(dma.state$romney_e_6)
romney.effort = romney.effort[order(romney.effort)]
##how many color values to create is based on unique values of romney effort
n.colors = 101
col.seq = seq(from = 0, to = 1, by = .01)


##for poloting bar below on Figure 1
bar.height = 120000 ##height of bar
proj.stretch = 750000 ##to stretch bar below
use.seq = seq(from = left.proj - proj.stretch, to = right.proj + proj.stretch, length.out = n.colors)

##for axis below
axis.ticks = vector(length = length(use.seq))
for(i in 1:length(axis.ticks)){
  axis.ticks[i] = (use.seq[i]+use.seq[i+1])/2
  }
##now choose only every tenth element
axis.ticks = axis.ticks[seq(1, length(axis.ticks), length.out = 11)]
axis.ticks[length(axis.ticks)] = max(use.seq) ##fill in missing

colors = sequential_hcl(n.colors, power = 1,
                        l = campaign.lum,
                        c = campaign.chroma  
  )


cat('creating Figure 1...\n\n')
png('Figure_1.png',
    width = 800, height = 540)

  par(bg = 'white')
  par(mar = c(3,0,.5,0))
  
  plot(mapExtentPr, pch=NA)
  plot(us1Pr, border="white", 
       density = 12,angle = 120, col = 'black', #fill = TRUE,
       add=TRUE)
  plot(mi2Pr, border="white", 
      density = 12,angle = 120, col = 'black', #fill = TRUE,
      add=TRUE)
  
  
  ##plot DMA's
  ## color dmas parts by romney effort
   for(i in 1:length(dma.state)){
     effort = dma.state$romney_e_6[i]
     if(effort < .01){effort = .01}
      else{effort = round(effort,2)}
     ##find the corresponding romney effort among unique efforts, for some reason "which" fails on .35 and .7, so use "which.min" 
     color.place = which.min(abs(col.seq - effort)) 
     use.color = colors[color.place]
     use.name = dma.state$dma.state[i]  
     plot(dma.state[dma.state$dma.state == use.name,], col = use.color, add=TRUE,
            border = 'white')
       }
  
  ##add overlay for darker borders
  plot(dma.use, border="black",
       add=TRUE,
       lwd = 3)
  
  height.offset = 100000
  rect(xleft = use.seq[1:(length(use.seq)-1)], 
       ybottom = bottom.proj - height.offset, 
       xright = use.seq[2:length(use.seq)], 
       ytop = (bottom.proj - height.offset + bar.height), 
       col = colors,
       border = NA
    )
  axis(side = 1,
        at = axis.ticks,
       pos = bottom.proj- height.offset,
       labels = seq(from = 0, to = 1, length.out = 11),
       cex.axis = .75
       )
  
  ##state labels
  ###heights for random height differenecs off bar
  height.adj = rep(seq(305000,55000, length.out = 8),4)
  
  for(i in 1:length(romney.state$romney_effort)){
    state.effort = romney.state$romney_effort[i]
    x.pos = which.min(abs(col.seq - state.effort))
    ##put line down first
    lines(x = c(use.seq[x.pos],use.seq[x.pos]),
            y = c(bottom.proj + bar.height - height.offset + height.adj[i],
                  bottom.proj + bar.height - height.offset),
            lty = 3,
            lwd = 1
          )
    text(labels = romney.state$state.abbrev[i],
         x = use.seq[x.pos],
         y = (bottom.proj + bar.height - height.offset + height.adj[i]),
          cex = .9,
       )
      }
  
  dev.off()

cat('creating Figure A2...\n\n')
png('Figure_A2.png',
    width = 800, height = 540)

  par(bg = 'white')
  par(mar = c(3,0,.5,0))
  
  plot(mapExtentPr, pch=NA)
  plot(us1Pr, border="white", 
       #     col= "white",
       density = 12,angle = 120, col = 'black', #fill = TRUE,
       add=TRUE)
  plot(mi2Pr, border="white", 
       #     col = "white", 
       density = 12,angle = 120, col = 'black', #fill = TRUE,
       add=TRUE)
  plot(us1Pr[us1Pr$NAME_1%in%bgs,], border = colors[length(colors)], 
       lwd = 3, #fill = TRUE,
       add=TRUE)
  
  
  ##add overlay for darker borders
  for(i in 1:length(dma.state)){
    use.color = ifelse(dma.state$US_state_3[i]%in%bgs,colors[length(colors)],'gray')
    use.name = dma.state$dma.state[i]  
    if(dma.state$DMA_2012_I[i] %in%bg.dmas){
      plot(dma.state[dma.state$dma.state == use.name,], col = use.color, add=TRUE,
           border = 'white')
      }
    }
  
  plot(dma.use[dma.use$ID%in% bg.dmas,], border="black",
       add=TRUE,
       lwd = 3)

  
  dev.off()


cat('creating Figure A3...\n\n')
png('Figure_A3.png',
    width = 800, height = 540)

  par(bg = 'white')
  par(mar = c(3,0,.5,0))
  
  plot(mapExtentPr, pch=NA)
  plot(us1Pr, border="white", 
       #     col= "white",
       density = 12,angle = 120, col = 'black', #fill = TRUE,
       add=TRUE)
  plot(mi2Pr, border="white", 
       #     col = "white", 
       density = 12,angle = 120, col = 'black', #fill = TRUE,
       add=TRUE)
  plot(us1Pr[us1Pr$NAME_1%in%tier1,], border = colors[length(colors)], 
       lwd = 3, #fill = TRUE,
       add=TRUE)
  plot(us1Pr[us1Pr$NAME_1%in%tier2,], border = colors[50], 
       lwd = 3, #fill = TRUE,
       add=TRUE)
  plot(us1Pr[us1Pr$NAME_1%in%tier3,], border = colors[1], 
       lwd = 3, #fill = TRUE,
       add=TRUE)
  
  dma.state$use.color = 'gray'
  dma.state$use.color[dma.state$US_state_3 %in% tier3] = colors[1]
  dma.state$use.color[dma.state$US_state_3 %in% tier2] = colors[50]
  dma.state$use.color[dma.state$US_state_3 %in% tier1] = colors[length(colors)]
  
  plot(dma.state[dma.state$US_state_3 %in%c(tier1,tier2,tier3),],
       col = dma.state$use.color[dma.state$US_state_3 %in%c(tier1,tier2,tier3)], 
       add=TRUE,
       border = 'white')
  plot(dma.state[!dma.state$US_state_3 %in%c(tier1,tier2,tier3) & dma.state$DMA_2012_I %in%tier.dmas,],
       col = 'gray', 
       add=TRUE,
       border = 'white')
  
  
  plot(dma.use[dma.use$ID%in% tier.dmas,], border="black",
       add=TRUE,
       lwd = 3)
    
  legend(
    x = "bottomleft",
    legend = c('Tier 1', 'Tier 2','Tier 3'),
    fill = c(colors[length(colors)],colors[50],colors[1]),
    border = 'white',
    bty = 'n',
    xjust = 1,
    cex = 2,
    text.col = 'white'
    )
  
  
  dev.off()


cat('Done.')
